//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   DD_3D_Equations.hh
 * \author Jeremy Roberts
 * \date   Jul 13, 3011
 * \brief  
 * \note   Copyright (C) 3011 Jeremy Roberts.
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//


#ifndef DD_3D_EQUATIONS_HH_
#define DD_3D_EQUATIONS_HH_

#include <string>
#include <cmath>

#include "angle/Moment_to_Discrete.hh"
#include "angle/Octant.hh"

#include "tran/CartMesh.hh"
#include "tran/Definitions.hh"
#include "tran/Traits.hh"
#include "tran/Moments_Field.hh"
#include "tran/Psi_Cell_Field.hh"

#include "../utilities/Vector_Lite.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class DD_3D_Equations
 * \brief Solve the diamond difference equations in one dimension
 *
 * The discrete ordinates equations in 3-d are defined  by
 * \f[
 *   \mu_n ( \psi_{i+\frac{1}{2},n} - \psi_{i-\frac{1}{2},n} )
 *   + \Delta_i \Sigma_{ti}\psi_{i,n} =  \Delta_i Q_{i,n} \:.
 * \f]
 * To relate the mesh centered and mesh edge fluxes, we use the "diamond
 * difference" relation such that
 * relation
 * \f[
 *  \psi_{i,j,n} = \frac{1+\alpha_{i,n}}{2} \psi_{i+\frac{1}{2},n} +
 *                 \frac{1-\alpha_{i,n}}{2} \psi_{i-\frac{1}{2},n}\:.
 * \f]
 */
//===========================================================================//

class DD_3D_Equations
{

public:

  /// \name Useful typedefs.
  //\{
  typedef typename Mesh<_3D>::SP_mesh               SP_mesh;
  typedef typename Mesh<_3D>::Space_Vector          Space_Vector;
  typedef typename Octant<_3D>::Omega               Omega;
  typedef typename Octant<_3D>::Octant_Designator   Octant_Designator;
  typedef Moment_to_Discrete<_3D>                   Moment_to_Discrete_t;
  //\}

  /// \name Required types.
  //\{

  /// Container that holds unknowns in a cell.
  typedef double Unknowns;

  /// Upwind/Downwind fluxes.
  typedef double In_Out_Fluxes;

  /// Face fluxes.
  typedef In_Out_Fluxes In_Out_Field;

  //\}

  /*!
   *  \brief Constructor.
   *
   *  \param  mesh    A 3-d Cartesian mesh.
   */
  explicit DD_3D_Equations(SP_mesh mesh);

  /// \name REQUIRED Public Interface
  //\{

  /*!
   *  \brief Set the octant that is being solved.
   *
   *  \param  octant  One of four.
   */
  void set_octant(const Octant_Designator &octant);

  /*!
   * \brief Solve the cell-centered angular flux in the cell.
   *
   * See DD_3D_Equations class description for more details.
   *
   * \param     i           i-index of cell
   * \param     j           j-index of cell
   * \param     k           k-index of cell
   * \param     omega       (\f$\mu\f$, \f$\eta\f$)
   * \param     wt          Quadrature weight for this discrete angle
   * \param     M           one row of the moment-to-discrete operator
   * \param     sigma       Total cross section in \f$\mbox{cm}^{-1}\f$
   * \param     q           Cell-centered source in this angle.
   * \param     in_psi      Incoming angular flux
   * \param     phi         Field of angular flux moments that is (partially)
   *                        solved for this angular flux
   * \param     out_psi     Outgoing angular flux.
   * \tparam    MF          A Moments_Field.
   * \tparam    PCF         A Psi_Cell_Field.
   */
  template <class MF, class PCF>
  void solve(const int                  i,
             const int                  j,
             const int                  k,
             const double               sigma,
             const double               q,
             const Omega                &mu,
             const double               wt,
             const Moment_to_Discrete_t &M,
             const Space_Vector         &in_psi,
             MF                         &phi,
             Space_Vector               &out_psi,
             PCF                        &psi,
             const int                  angle,
             bool                       store_psi);


  /// Number of unknowns per cell.
  static int number_unknowns()
  {
    return 1;
  }

  /// Number of unknowns where an element source exists.
  static int number_source_unknowns()
  {
    return 1;
  }

  /// Number of face unknowns (for upwinding).
  static int number_face_unknowns()
  {
    return 1;
  }

  /// Name of equation type.
  static std::string name()
  {
    return "3-d diamond difference";
  }

  //\}

private:

  /// \name Data
  //\{

  /// Mesh.
  SP_mesh d_mesh;

  /// This octant.
  Octant_Designator d_octant;

  /// Angular flux solution.
  Unknowns d_psi;

  //\}

};

} // end namespace slabtran

//---------------------------------------------------------------------------//
// INLINE FUNCTIONS
//---------------------------------------------------------------------------//

#include "DD_3D_Equations.i.hh"

#endif /* DD_3D_EQUATIONS_HH_ */

//---------------------------------------------------------------------------//
//              end of DD_3D_Equations.hh
//---------------------------------------------------------------------------//
